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ABSTRACT 


This  report  forms  the  user's  guide  for  Version  1.1  of  SOL/NPSOL,  a  set  of  Fortran  subroutines 
designed  to  minimise  an  arbitrary  smooth  function  subject  to  constraints,  which  may  include 
ample  bounds  on  the  variables,  linear  constraints  and  smooth  nonlinear  constraints.  (NPSOL 
may  also  be  used  for  unconstrained,  bound*constrained  and  linearly  constrained  optimisation.) 
The  user  must  provide  subroutines  that  define  the  objective  and  constraint  functions  and  their 
gradients.  All  matrices  are  treated  as  dense,  and  hence  NPSOL  is  not  intended  Tor  large  sparse 
problems. 

NPSOL  uses  a  sequential  quadratic  programming  (SQP)  algorithm,  in  which  the  search  direc¬ 
tion  is  the  solution  of  a  quadratic  programming  (QP)  subprobiem.  The  algorithm  treats  bounds, 
linear  constraints  and  nonlinear  constraints  separately.  The  Hessian  of  each  QP  subproblem 
is  a  positive-definite  quasi-Newton  approximation  to  the  Hessian  of  an  augmented  Lagrangian 
function.  The  steplength  at  each  iteration  is  required  to  produce  a  sufficient  decrease  in  an  aug¬ 
mented  Lagrangian  merit  function.  Each  QP  subproblem  is  solved  using  a  quadratic  programming 
package  with  several  features  that  improve  the  efficiency  of  an  SQP  algorithm.  ~ _ 


*The  package  SOL/NPSOL  is  available  from  the  Office  of  Technology  Licensing,  105  Encina  Hall, 
Stanford  University,  Stanford,  California,  94305. 
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1.  PURPOSE 


SOL/NPSOL  is  a  collection  of  Fortran  subroutines  designed  tc  solve  the  nonlinear  program¬ 
ming  problem  —  the  minimization  of  a  smooth  nonlinear  function  subject  to  a  set  of  constraints 
on  the  variables.  The  problem  is  assumed  to  be  stated  in  the  following  form: 


NP 

minimize 

*■(*) 

*eR" 

X 

| 

subject  to 

l  <  ' 

Alx 

c(x) 

<  «, 

where  F[x)  is  a  smooth  nonlinear  function,  AL  is  a  constant  matrix,  and  c(z)  is  a  vector  of  smooth 
nonlinear  constraint  functions.  The  matrix  AL  and  the  vector  c(z)  may  be  empty.  Note  that  upper 
and  lower  bounds  are  specified  for  all  the  variables  and  for  all  the  constraints.  This  form  allows 
full  generality  in  specifying  other  types  of  constraints.  In  particular,  the  t-th  constraint  may 
be  defined  as  an  equality  by  setting  U  =  ut.  If  certain  bounds  are  not  present,  the  associated 
elements  of  t  or  u  can  be  set  to  special  values  that  will  be  treated  as  —  oo  or  +oo. 

If  no  nonlinear  constraints  ?re  present,  it  is  generally  more  efficient  to  use  a  package  specifically 
designed  for  linearly  constrained  problems.  In  particular,  when  F  is  linear  or  quadratic,  the 
LPSOL  or  QPSOL  packages  should  be  used  (Gill  et  a/.,  1983a);  for  a  general  function  F  with 
only  linear  constraints,  the  LCSOL  package  is  appropriate  (Gill  et  a/.,  1983c).  If  the  problem 
is  large  and  sparse,  the  MINOS/ AUGMENTED  package  (Murtagh  and  Saunders,  1980,  1982) 
should  be  used,  since  NPSOL  treats  all  matrices  as  dense. 

The  user  must  supply  an  initial  estimate  of  the  solution  to  NP,  and  subroutines  that  define 
F{x),  c(x)  and  their  first  derivatives.  The  level  of  printed  output  is  controlled  by  the  user  (see 
the  parameter  MSGLVL  in  Section  4). 

NPSOL  is  based  on  subroutines  from  Version  3.1  of  the  SOL/QPSOL  quadratic  programming 
package;  the  documentation  of  this  version  of  QPSOL  (Gill  et  a/.,  1983a)  should  be  consulted  in 
conjunction  with  this  report.  NPSOL  contains  approximately  9000  lines  of  ANSI  (1966)  Standard 
Fortran,  of  which  47%  are  comments. 
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The  method  used  to  solve  NP  is  a  sequential  quadratic  programming  (SQP)  method.  SQP 
methods  were  popularized  mainly  by  Biggs  (1972),  Han  (1976)  and  Powell  (1977);  for  an  overview, 
see,  e.g.,  Fletcher  (1981),  Gill,  Murray  and  Wright  (1981)  and  Powell  (1982).  Let  xq  denote  the 
initial  estimate  of  the  solution.  During  the  fc-th  “major  iteration”  of  NPSOL  (k  =  0,1,...),  a 
new  estimate  is  defined  by 

**+ 1  =  *fc  +  ttfcPfei 

where  the  vector  pk  is  the  solution  of  a  QP  subproblem,  to  be  described  below.  The  positive 
scalar  a*  is  chosen  to  produce  a  sufficient  decrease  in  an  augmented  Lagrangian  merit  function 
(see  Schittkowski,  1981);  the  procedure  that  determines  a*  is  called  the  line  search. 

The  QP  subproblem  that  defines  p*  is  of  the  form 


QP 


minimise 

pea- 


subject  to 


The  vector  g  in  QP  is  the  gradient  of  F  at  z*.  The  matrix  H  is  a  positive-definite  quasi-Newton 
approximation  to  the  Hessian  of  an  augmented  Lagrangian  function.  It  is  represented  as  H  = 
RtR,  where  R  is  upper  triangular,  and  is  updated  after  every  major  iteration. 

Let  mL  denote  the  number  of  linear  constraints  (the  number  of  rows  in  AL),  and  let  mN 
denote  the  number  of  nonlinear  constraints  (the  dimension  of  c(x)).  The  matrix  A  in  QP  has 
mL  +  mN  rows,  and  is  defined  as 


A  = 


where  AH  is  the  Jacobian  matrix  of  c(z)  evaluated  at  z*.  Let  t  in  NP  be  partitioned  into  three 
sections:  the  first  n  components  (denoted  by  tB),  corresponding  to  the  bound  constraints;  the 
next  mL  components  (denoted  by  tL),  corresponding  to  the  linear  constraints;  and  the  last  mN 
components  (denoted  by  lN),  corresponding  to  the  nonlinear  constraints.  The  vector  i  in  QP  is 
partitioned  in  the  same  way,  and  is  defined  as 


tg  —  f-B  *fc»  "f-L  —  and  In  —  tfi  Cfc, 


where  c*  is  c(z)  evaluated  at  z^.  The  vector  u  is  defined  in  an  analogous  fashion. 

In  general,  solving  the  subproblcm  QP  for  pk  is  itself  an  iterative  procedure.  Hence,  a 
“minor  iteration”  of  NPSOL  corresponds  to  an  iteration  within  the  QP  algorithm.  Note  that 
the  functions  Fix)  and  c(z)  are  not  evaluated  during  the  solution  of  the  subproblem.  The  total 
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number  of  function  evaluations  required  to  solve  a  well-behaved  problem  will  usually  be  similar 
to  the  number  of  major  iterations. 

The  problem  QP  is  solved  using  subroutines  from  the  SOL/QPSOL  package,  which  is  described 
in  detail  in  Gill  et  a /.  (1983a),  and  was  specifically  designed  to  be  used  within  an  SQP  algorithm 
for  nonlinear  programming.  In  particular,  two  common  difficulties  associated  with  SQP  methods 
are  alleviated  by  certain  features  of  the  QPSOL  subroutines. 

First,  it  may  happen  that  the  QP  subproblem  is  infeasible,  yet  feasible  points  exist  with 
respect  to  the  nonlinear  constraints.  (Throughout  this  report,  we  assume  that  “feasibility”  is 
defined  by  a  set  of  tolerances  provided  by  the  user  in  the  array  FEATOL;  see  Section  4.)  The 
strategy  used  by  NPSOL  to  treat  an  infeasible  subproblem  is  the  following.  If  there  is  no  feasible 
point  with  respect  to  the  bounds  and  linear  constraints  of  the  original  problem,  the  infeasibility  is 
inherent  in  the  problem,  and  hence  NPSOL  terminates.  Otherwise,  the  infeasibility  results  from 
the  linearized  nonlinear  constraints;  the  least  infeasible  point  is  then  computed,  the  appropriate 
constraint  bounds  are  (temporarily)  relaxed,  and  a  relaxed  quadratic  program  is  solved  for  p*. 

Second,  it  is  useful  in  an  SQP  algorithm  to  be  able  to  use  the  prediction  of  the  active  set 
from  each  QP  subproblem  to  solve  the  next  subproblem  more  efficiently.  This  benefit  is  achieved 
in  NPSOL  by  a  “hot  start”  feature  that  allows  the  initial  working  set  and  part  of  its  factorization 
to  be  specified.  Within  NPSOL,  the  prediction  of  the  active  set  from  one  QP  subproblem  is  used 
as  the  “hot  start”  estimate  of  the  working  set  for  the  next  QP.  In  practice,  this  means  that  the 
QP  subproblems  near  the  solution  reach  optimality  in  only  one  iteration.  Furthermore,  separate 
treatment  of  linear  constraints  means  that  it  is  usually  possible  to  save  work  in  performing  the 
factorization  of  the  working  set  at  the  beginning  of  the  QP  (since  the  rows  of  A  corresponding  to 
the  linear  constraints  are  unchanged). 

The  algorithm  used  in  NPSOL  will  be  discussed  in  a  forthcoming  report.  Details  of  the 
algorithm  of  QPSOL  are  given  in  Gill  et  at.  (1983b). 
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SUBROUTINE  NPSOL  (  UMAX,  MSGLVL,  N, 

NCLIN,  NCNLN ,  NCTOTL,  NROVA,  NROWJ,  NROWR, 
BIGBND,  EPSAF,  ETA,  FTOL, 

A,  BL.  BU,  FEATOL, 

CONFUN,  OBJFUN,  COLD.  FEALIN,  ORTHOG, 
INFORM,  ITER,  ISTATE, 

C,  CJAC,  CLAUD A,  OBJF,  OBJGRD,  R,  X, 

IV,  LENIN,  W.  LENV  ) 


EXTERNAL 

CONFUN,  OBJFUN 

LOGICAL 

COLD.  FEALIN,  ORTHOG 

INTEGER 

UMAX,  MSGLVL,  N.  NCLIN.  NCNLN,  NCTOTL. 

NROffA,  NROWJ,  NROWR,  INFORM.  ITER.  LENIN,  LENW 

INTEGER 

ISTATE (NCTOTL) .  IN (LENIN) 

REAL 

BIGBND.  EPSAF.  ETA.  FTOL.  OBJF 

REAL 

A(NROWA.N),  BL (NCTOTL) ,  BU (NCTOTL) ,  FEATOL (NCTOTL) 
C (NROWJ) ,  CJAC (NROWJ, N),  CLAMDA (NCTOTL) , 

OBJGRD (N) ,  R (NROWR, N) ,  X(N),  N(LENN) 

Note:  Here  and  elsewhere,  the  specification  of  a  parameter  as  REAL  should  be  interpreted  as 
working  precision,  which  may  be  DOUBLE  PRECISION  in  some  circumstances. 
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ITMAX  is  an  upper  bound  on  the  number  of  major  iterations  to  be  performed.  Unless  the 
problem  is  known  to  be  exceptionally  difficult,  a  sensible  initial  choice  for  ITMAX  is  50. 

MSGLVL  indicates  the  amount  of  intermediate  output  desired  (see  Section  9  for  a  description  of 
the  printout).  All  output  is  written  to  the  file  number  NOUT  (see  subroutine  HCHPAR  in 
Section  11).  MSGLVL  is  interpreted  as  a  four-digit  number.  Its  first  two  digits  indicate 
the  level  of  intermediate  output  from  the  quadratic  programming  routines;  the  second 
two  digits  indicate  the  level  of  intermediate  output  from  NPSOL.  The  QP  printout 
levels  are  defined  in  Gill  et  a 1.  (1983a);  if  MSGLVL  <  100,  there  is  no  QP  output.  When 
the  last  two  digits  of  MSGLVL  >  10,  each  level  includes  the  printout  from  all  lower  levels. 
The  printout  corresponding  to  each  value  of  the  last  two  digits  of  MSGLVL  is  defined  aa 
follows: 

Value  Definition 

0  No  output. 

1  The  final  solution  only. 

5  One  brief  line  of  output  for  each  major  iteration  (no  printout  of  the 
final  solution). 

>  10  The  final  solution  and  one  brief  line  of  output  for  each  major  iteration. 

>  15  At  each  iteration,  the  arrays  X  and  I  STATE,  and  the  indices  of  the  free 

variables. 

>  20  At  each  iteration,  the  nonlinear  constraint  values  (the  array  C),  the 

linear  constraint  values  (ALx),  and  estimates  of  the  Lagrange  multi¬ 
pliers. 

>  30  At  each  iteration,  the  diagonal  elements  of  the  matrix  T  associated 

with  the  TQ  factorisation  of  the  working  set,  and  the  diagonals  of  the 
matrix  R  (the  Cholesky  factor  of  the  Hessian  approximation). 

>  80  Debug  output  from  NPSOL. 

99  Debug  output  from  the  line  search. 

For  example,  MSGLVL  =  10  will  produce  a  summary  of  results  for  each  major  iteration 
and  a  full  printout  of  the  final  solution;  MPCLVL  =  510  will  produce  the  same  printout, 
as  well  as  a  summary  of  cac»  minor  (O'  .tcration. 

R  is  the  number  of  variables,  i.e.,  the  dimension  of  X  (N  must  be  positive). 
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NCLIN  is  the  number  of  general  linear  constraints  in  the  problem  (NCLIN  may  be  zero). 

NCNLN  is  the  number  of  nonlinear  constraints  in  the  problem  (NCNLN  may  be  zero). 

NCTOTL  must  be  set  to  N  +  NCLIN  +  NCNLN. 

NROWA  is  the  declared  row  dimension  of  the  array  A  (NROWA  must  be  at  least  1  and  at  least 
NCLIN). 

NROWJ  is  the  declared  row  dimension  of  the  array  CJAC  and  the  length  of  the  array  C  (NROWJ 
must  be  at  least  1  and  at  least  NCNLN). 

NROWR  is  the  declared  row  dimension  of  the  array  R  (NROWR  must  be  at  least  N). 

BIGBND  is  a  positive  real  variable  whose  magnitude  denotes  an  “infinite”  component  of  l  and 
v.  Any  upper  bound  greater  than  or  equal  to  BIGBND  will  be  regarded  as  plus  infinity 
(and  similarly  for  a  lower  bound  less  than  or  equal  to  —BIGBND). 

EPSAF  is  a  positive  quantity  that  should  be  a  good  bound  on  the  absolute  error  in  computing 
F(x)  at  the  initial  point.  For  many  simple  functions,  EPSAF  is  of  the  order  of  cM|/'’(z)|, 
where  is  the  machine  precision.  A  discussion  of  EPSAF  is  given  in  Chapter  8  of  Gill, 
Murray  and  Wright  (1981). 

ETA  is  a  number  satisfying  0  <  ETA  <  1,  which  controls  how  accurately  the  value  a* 
approximates  a  univariate  minimum  of  the  merit  function  along  p*  (the  smaller  the 
value  of  ETA,  the  more  accurate  the  line  search).  The  recommended  value  of  ETA  for 
nonlinearly  constrained  problems  is  0.9,  which  corresponds  to  a  relaxed  line  search. 
If  the  problem  is  unconstrained,  bound-constrained,  or  linearly  constrained,  a  smaller 
value  of  ETA  will  tend  to  require  more  function  evaluations,  but  fewer  major  iterations. 

FTOL  is  a  positive  tolerance  (FTOL  <  1)  that  indicates  the  number  of  figures  of  accuracy 
desired  in  the  objective  function  at  the  solution.  For  example,  if  FTOL  is  10-6  and 
NPSOL  terminates  successfully,  the  computed  solution  should  have  approximately  six 
correct  figures  in  F.  FTOL  should  never  be  less  than  machine  precision. 

A  is  a  real  array  of  declared  dimension  (NROWA,  N),  corresponding  to  A,,  in  the  problem 

specification  NP  (Section  1).  The  i-th  row  of  A,  »  =  1  to  NCLIN,  contains  the  coefficients 
of  the  i-th  general  linear  constraint.  If  NCLIN  is  zero,  A  is  not  accessed. 
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BL  is  a  real  array  of  dimension  NCTOTL  that  contains  the  lower  bounds  for  all  the  constraints, 

in  the  following  order  (which  is  also  observed  for  BU,  CLAMDA,  FEATOL  and  ISTATE).  The 
first  N  elements  of  BL  contain  the  lower  bounds  on  the  variables.  If  NCLIN  >  0,  the  next 
NCLIN  elements  of  BL  contain  the  lower  bounds  for  the  general  linear  constraints.  If 
NCNLN  >  0,  the  next  NCNLN  elements  of  BL  contain  the  lower  bounds  for  the  nonlinear 
constraints.  In  order  for  the  problem  specification  to  be  meaningful,  it  is  required  that 
BL(»  <  BU(i)  for  all  j.  To  specify  a  non-existent  lower  bound  for  the  j- th  constraint 
(i.e.,  lj  =  — oo),  the  value  used  must  satisfy  BL(j)  <  —  BIGBND.  To  specify  the  j-th 
constraint  as  an  equality ,  the  user  must  set  BL(j')  =  BU(j)  =  /J,  say  where  |/3(  < 
BIGBND. 

BU  is  a  real  array  of  dimension  NCTOTL  that  contains  the  upper  bounds  for  all  the  con¬ 

straints,  in  the  same  order  described  above  for  BL.  To  specify  a  non-existent  upper 
bound  (i.e.,  uy  =  +00),  the  value  used  must  satisfy  BU(j)  >  BIGBND. 

FEATOL  is  a  real  array  of  dimension  NCTOTL  containing  positive  tolerances  that  define  the 
maximum  permissible  violation  in  each  constraint  in  order  for  a  point  to  be  considered 
feasible,  i.e.  constraint  j  is  considered  satisfied  if  its  violation  does  not  exceed  FEATOL(j). 
The  ordering  of  the  components  of  FEATOL  is  the  same  as  that  described  above  under  BL. 
Note  that  FEATOL(j)  is  a  bound  on  the  absolute  acceptable  violation.  For  example,  if  the 
data  defining  the  constraints  are  of  order  unity  and  are  correct  to  about  6  decimal  digits, 
it  would  be  appropriate  to  choose  FEATOL(j)  as  10~8  for  all  relevant  j.  In  general,  the 
elements  of  FEATOL  should  be  chosen  as  the  largest  possible  acceptable  values,  since  the 
algorithm  of  NPSOL  becomes  less  likely  to  encounter  difficulties  with  ill-conditioning 
and  degeneracy  as  the  components  of  FEATOL  increase.  A  warning  message  is  printed 
if  any  component  of  FEATOL  is  less  than  machine  precision;  the  user  must  not  set  any 
element  of  FEATOL  to  zero.  A  detailed  discussion  of  FEATOL  is  given  in  Gill  et  al.  (1983b). 

CONFUN  is  the  name  of  a  subroutine  that  calculates  the  vector  c(x )  of  nonlinear  constraint 
functions  and  its  Jacobian  for  a  specified  n-vector  x.  CONFUN  must  be  declared  as 
EXTERNAL  in  the  routine  that  calls  NPSOL.  IT  there  are  no  nonlinear  constraints  (NCNLN 
=  0),  CONFUN  will  never  be  called  by  NPSOL.  IT  there  arc  nonlinear  constraints,  NPSOL 
always  calls  CONFUN  and  OBJFUN  together,  in  that  order. 

The  specification  of  CONFUN  is: 

SUBROUTINE  CONFUN (  MODE,  NCNLN,  N,  NROWJ,  X.  C,  CJAC,  NSTATE  ) 

INTEGER  MODE,  NCNLN,  N,  NROWJ,  NSTATE 

REAL  X(N) ,  C (NROWJ) ,  CJAC (NROWJ, N) . 

The  actual  parameters  NCNLN,  N,  and  NROWJ  input  to  CONFUN  will  always  be  the 

same  Fortran  variables  as  those  input  to  NPSOL.  They  must  not  be  altered  by  CONFUN. 
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4.  INPUT  PARAMETERS 


MODE  is  a  flag  that  the  user  may  set  within  CONFUN  to  indicate  a  failure  in  the 
evaluation  of  the  nonlinear  constraints.  On  entry  to  CONFUN,  MODE  is  always  nonnegativc. 
If  MODE  is  negative  on  exit  from  CONFUN,  the  execution  of  NPSOL  will  be  terminated  with 
INFORM  set  to  MODE. 

X  contains  the  vector  of  variables  z  at  which  the  constraint  functions  are  to  be 
evaluated.  The  elements  of  X  must  not  be  altered  by  CONFUN. 

C  should  contain  the  nonlinear  constraint  values  ct(z),  t  =  1  to  NCNLN,  on  exit  from 
CONFUN. 

CJAC  should  contain  the  Jacobian  matrix  of  the  nonlinear  constraint  functions  on 
exit  from  CONFUN.  The  i-th  row  of  CJAC  contains  the  gradient  of  the  i-th  nonlinear 
constraint,  i.e.  CJAC(i,jf)  is  the  partial  derivative  of  c,  with  respect  to  x},  i  =  1  to 
NCNLN,  j  =  1  to  N.  If  CJAC  contains  any  constant  elements,  a  saving  in  computation 
can  be  made  by  setting  them  one  time  only,  when  NSTATE  =  1  (see  below). 

NSTATE  is  set  to  one  by  NPSOL  on  the  first  call  of  CONFUN,  and  is  zero  for  all 
subsequent  calls.  Thus,  if  the  user  wishes,  NSTATE  may  be  tested  within  CONFUN  in 
order  to  perform  certain  calculations  one  time  only.  For  example,  the  user  may  read 
data  or  initialize  COMMON  blocks  when  NSTATE  =  1.  In  addition,  the  constant  elements 
of  CJAC  can  be  set  in  CONFUN  when  NSTATE  =  1,  and  need  not  be  defined  on  subsequent 
calls. 

OBJFUN  is  the  name  of  a  subroutine  that  calculates  the  objective  function  F(x)  and  its  gradient 
for  a  specified  n-vector  z.  OBJFUN  must  be  declared  as  EXTERNAL  in  the  routine  that 
NPSOL. 

The  specification  of  OBJFUN  is: 

SUBROUTINE  OBJFUN (  MODE.  N.  X,  OBJF,  OBJGRD.  NSTATE  ) 


INTEGER 


MODE,  N,  NSTATE 
OBJF.  X(N) ,  OBJGRD (N) . 


The  actual  parameter  N  input  to  OBJFUN  will  always  be  the  same  Fortran  variable 
as  that  input  to  NPSOL,  and  must  not  be  altered  by  OBJFUN. 

MODE  is  a  flag  that  the  user  may  set  within  OBJFUN  to  indicate  a  failure  in  the 
evaluation  of  the  objective  function.  On  entry  to  OBJFUN,  MODE  is  always  nonnegativc. 
If  MODE  is  negative  on  exit  from  OBJFUN,  the  execution  of  NPSOL  is  terminated  with 
INFORM  set  to  MODE. 

X  contains  the  vector  of  variables  z  at  which  the  objective  function  is  to  be  evaluated. 
The  X  array  must  not  be  altered  by  OBJFUN. 

OBJF  should  contain  the  value  of  the  objective  function  /'’(z)  on  exit  from  OBJFUN. 

OBJGRD  should  contain  the  gradient  vector  of  the  objective  function.  The  j-th 
component  of  OBJGRD  contains  the  partial  derivative  of  F  with  respect  to  the  j-lh 
variable. 


4.  INPUT  PARAMETERS 
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NSTATE  is  set  to  one  by  NPSOL  on  the  first  call  of  OBJFUN,  and  to  zero  on  ail 
subsequent  calls.  Thus,  if  the  user  wishes,  NSTATE  may  be  tested  in  order  to  perform 
certain  calculations  only  on  the  first  call  of  OBJFUN  —  e.g.,  read  data  or  initialize  COMMON 
blocks.  Note  that  if  there  are  any  nonlinear  constraints,  CONFUN  and  OBJFUN  are  always 
called  together,  in  that  order. 

COLD  is  a  logical  variable  that  indicates  whether  the  user  has  specified  an  initial  estimate  of 
the  active  set  of  constraints.  If  COLD  is  .TRUE.,  the  initial  working  set  is  determined  by 
the  first  QP  subproblem.  If  COLD  is  .FALSE,  (a  “warm  start”),  the  user  must  define  the 
array  ISTATE  (which  gives  the  status  of  each  constraint  with  respect  to  the  working  set) 
and  the  matrix  R  (the  Cholesky  factor  of  the  initial  Hessian  approximation).  The  warm 
start  option  is  particularly  useful  when  NPSOL  is  restarted  at  the  point  where  an  earlier 
run  terminated. 

FEALIN  is  a  logical  variable  that  indicates  whether  the  starting  point  for  the  SQP  method 
should  first  be  made  feasible  with  respect  to  the  bounds  and  linear  constraints  of  NP. 
If  FEALIN  is  .TRUE.,  the  algorithm  will  determine  (if  possible)  a  point  that  is  feasible 
with  respect  to  the  bounds  and  linear  constraints  before  beginning  the  SQP  iterations 
(where  “feasible"  is  defined  by  the  array  FEATOL;  see  above).  This  setting  of  FEALIN 
ensures  that  a//  iterates  within  the  SQP  algorithm  will  be  feasible  with  respect  to  the 
bounds  and  linear  constraints  (this  may  be  essential  in  certain  applications).  If  FEALIN 
is  .FALSE.,  the  SQP  method  will  begin  with  the  user-specified  initial  value  of  X.  In  this 
case,  the  iterates  will  not  necessarily  be  feasible  with  respect  to  the  linear  constraints 
of  the  original  problem  (unless  the  original  point  is  feasible).  In  general,  we  recommend 
a  value  of  .TRUE,  for  FEALIN. 

ORTHOG  is  a  logical  variable  that  indicates  whether  orthogonal  transformations  will  be  used  in 
the  QP  algorithm  to  compute  and  update  the  TQ  factorization  of  the  working  set 

AQ  =  (0  T), 

where  A  is  a  submalrix  of  A  and  T  is  reverse-triangular  (see  (till  ct  a/.,  1982).  IT  ORTHOG 
is  .TRUE.,  the  TQ  factorization  is  computed  using  Householder  reflections  and  plane 
rotations,  and  the  matrix  Q  is  orthogonal.  If  ORTHOG  is  .FALSE.,  stabilized  elementary 
transformations  arc  used  to  maintain  the  factorization,  and  Q  is  not  orthogonal.  A  rule 
of  thumb  in  making  the  choice  is  that  orthogonal  transformations  require  more  work, 
but  provide  greater  numerical  stability.  Thus,  we  recommend  setting  ORTHOG  to  .TRUE, 
in  any  of  the  following  situations:  the  problem  is  reasonably  small;  the  functions  arc 
highly  nonlinear;  the  active  set  is  ill-conditioned;  or  the  time  required  to  compute  the 
TQ  factorization  is  not  significant  compared  to  the  evaluation  of  the  problem  functions. 


/  / 
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ISTATE  is  an  integer  array  of  dimension  NCTOTL  that  indicates  the  status  of  every  constraint 
with  respect  to  the  current  prediction  of  the  active  set.  The  ordering  of  ISTATE  is  the 
same  as  that  described  above  for  BL,  i.e.,  the  first  N  components  of  ISTATE  refer  to  the 
bounds  on  the  variables,  the  next  NCLIN  components  refer  to  the  linear  constraints,  and 
the  last  NCNLN  components  refer  to  the  nonlinear  constraints.  The  significance  of  each 
possible  value  of  ISTATE(j)  is  as  follows: 


ISTATE  (j) 
-2 

-1 

0 

1 

2 

3 


Meaning 

This  constraint  (or  its  linearization)  violates  its  lower  bound  by  more 
than  FEATOL(j)  in  a  QP  subproblem. 

This  constraint  (or  its  linearization)  violates  its  upper  bound  by  more 
than  FEATOL(j)  in  a  QP  subproblem. 

The  constraint  is  not  in  the  predicted  active  set. 

This  inequality  constraint  is  included  in  the  predicted  active  set  at  its 
lower  bound. 

This  inequality  constraint  is  included  in  the  predicted  active  set  at  its 
upper  bound. 

The  constraint  is  included  in  the  predicted  active  set  as  an  equality. 
This  value  of  ISTATE  can  occur  only  when  BL(j)  =  BU(j). 


If  COLD  =  .TRUE.,  ISTATE  need  not  be  set  by  the  user.  However,  when  COLD  is 
.FALSE.,  every  element  of  ISTATE  must  be  set  to  one  or  the  values  given  above  to  define 
a  suggested  prediction  of  the  active  set  (which  will  be  used  as  the  initial  working  set  in 
the  first  QP  subproblem).  The  most  likely  values  are: 


ISTATE(j) 


Meaning 


0  The  corresponding  constraint  should  not  be  in  the  initial  working  set. 

1  The  constraint  should  be  in  the  initial  working  set  at  its  lower  bound. 

2  The  constraint  should  be  in  the  initial  working  set  at  its  upper  bound. 

3  The  constraint  should  be  in  the  initial  working  set  as  an  equality.  This 
value  must  not  be  specified  unless  BL  (j)  =  BU(j').  The  input  values  1, 
2  or  3  of  ISTATE(j)  all  have  the  same  effect  when  BL(j)  =  B U(j). 


On  exit  from  NPSOL,  the  values  in  the  ISTATE  array  indicate  the  composition  of  the 
active  set  of  the  final  QP  subproblem. 
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R  is  a  real  array  of  declared  dimension  (NR0WR,N)  that  contains  the  upper-triangular 

Cholcsky  factor  of  the  current  approximation  of  the  Hessian  of  the  Lagrangian  function. 
If  COLD  is  .TRUE.,  the  array  R  need  not  be  initialized  by  the  user.  If  COLD  is  .FALSE.,  R 
must  contain  an  appropriate  upper-triangular  matrix. 

X  is  a  real  array  of  dimension  N  that  contains  the  current  estimate  of  the  solution.  On 

entry  to  NPSOL,  X  must  be  defined;  on  exit  from  NPSOL,  X  contains  the  final  estimate  of 
the  solution. 


6.  OUTPUT  PARAMETERS 
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8.  OUTPUT  PARAMETERS 

INFORM  is  an  integer  that  indicates  the  result  of  NPSOL.  (When  MSGLVL  >  0,  a  short  description 
of  INFORM  is  printed.)  The  possible  values  of  INFORM  are: 

INFORM  Definition 

<  0  The  user  has  set  MODE  to  this  negative  value  in  CONFUN  or  OBJFUN. 

0  X  satisfies  the  first-order  optimality  conditions,  i.e.,  the  projected  gra¬ 
dient  and  the  active  constraint  residuals  are  negligible,  and  the  La¬ 
grange  multipliers  indicate  optimality. 

1  No  feasible  point  could  be  found  for  the  linear  constraints  and  bounds. 

2  No  improved  point  for  the  merit  function  could  be  found  during  the 
final  line  search. 

3  The  limit  of  ITMAX  major  iterations  was  reached. 

4  Extremely  small  Lagrange  multipliers  could  not  be  resolved. 

5  A  descent  direction  for  the  merit  function  could  not  be  found. 

0  An  input  parameter  is  invalid. 

ITER  is  an  integer  that  gives  the  number  of  major  iterations  performed. 

C  is  a  real  array  of  dimension  NROVJ  that  contains  the  values  of  the  nonlinear  constraint 

functions  C(t),  *  =  1  to  NCNLN,  at  the  final  iterate.  If  NCNLN  =  0,  C  is  not  accessed  by 
NPSOL. 

CJAC  is  a  real  array  of  dimension  (NROVJ, N)  that  contains  the  Jacobian  matrix  of  the  nonlinear 

constraint  functions  at  the  final  iterate,  i.e.  CJAC(t,  j)  contains  the  partial  derivative  of 
the  t'-th  constraint  function  with  respect  to  the  j-th  variable,  i  =  1  to  NCNLN,  j  =  1 
to  N.  If  NCNLN  =  0,  CJAC  is  not  accessed  by  NPSOL.  (Sec  the  discussion  of  CJAC  under 
CONFUN  above.) 

CLAMDA  is  a  real  array  of  dimension  NCTOTL  that  contains  the  final  multiplier  estimate  for  every 
constraint  (i.e.,  the  multipliers  of  the  final  QP  subproblcm).  The  ordering  of  CLAMDA 
is  the  same  as  that  given  above  for  BL.  If  the  j-th  constraint  is  defined  as  “inactive" 
by  the  ISTATE  array,  CLAMDA(j)  should  be  *ero;  if  the  j-th  constraint  is  an  inequality 
active  at  its  lower  bound,  CLAMDA(j)  should  be  non-negative;  if  the  j-th  constraint  is  an 
inequality  active  at  its  upper  bound,  CLAMDA(j)  should  be  non-positive. 

OBJF  is  the  value  of  the  objective  function  F(z)  at  the  final  iterate. 

OBJGRD  is  a  real  array  of  dimension  N  that  contains  the  gradient  of  the  objective  function. 
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7.  WORKSPACE  PARAMETERS 


IV  is  an  integer  array  of  dimension  LENIV,  which  provides  integer  workspace  for  NPSOL. 
LENIW  is  the  dimension  of  IV,  and  must  be  at  least  2N. 

V  is  a  real  array  of  dimension  LENV,  which  provides  real  workspace  for  NPSOL. 

LENV  is  the  dimension  of  V,  and  must  be  at  least  2N2  +  N(NC0N  +  NROV J  +  6)  +  2NC0N  +  NROVA  + 

max(lON  +  2NC0N  +  NROVA  +  NROVJ,  5N  +  4NC0N),  where  NCON  =  max(l,  NCLIN  +  NCNLN). 
An  overestimate  of  this  number  is  2NS  +  N(NC0N  + NROVJ +  16)  +  6NCON+2NROVA+ NROVJ. 

If  MSGLVL  >  0,  the  amount  of  workspace  provided  and  the  amount  of  workspace  required 
are  printed.  As  an  alternative  to  computing  LENV  from  the  formula  given  above,  the  user  may 
prefer  to  obtain  an  appropriate  value  from  the  output  of  a  preliminary  run  with  a  positive  value 
of  MSGLVL  and  LENV  set  to  1  (NPSOL  will  then  terminate  with  INFORM  =  0). 


8.  AUXILIARY  SUBPROGRAMS  AND  LABELLED  COMMON 


SPSOL/IS 


«.  AUXILIARY  SUBPROGRAMS  AND  LABELLED  COMMON 


The  auxiliary  subroutines  used  by  NPSOL  may  be  divided  into  three  groups.  The  first  group 
includes  the  following  subroutines,  which  are  not  part  of  the  QP  package: 


GETPTC 


NPQPGN 
R1BFGS 


NPCORE 


NPRHO 

R1M0D. 


NPGETC 


NPSRCH 


NPGETF 
’R1 

NPTQ 


The  second  group  of  subroutines  —  those  used  by  the  QP  package  —  are: 


ADDCON 

ALLOC 

BDPERT 

BNDALF 

CHKDAT 

DELCON 

FINDP 

GETLAM 

LPBGST 

LPCORE 

LPCRSH 

LPDUMP 

LPGRAD 

LPPRT 

MO  VEX 

QPCHKP 

QPCOLR 

QPCORE 

QPCRSH 

QPDUMP 

QPGRAD 

QPPRT 

PRTSOL 

RSOLVE 

TQADD 

TSOLVE 

ZYPROD . 

NPSOL  also  uses  the  basic  linear  algebra  subroutines 


AXPY 

CONDVC 

COPYMX 

COPYVC 

DOT 

DSCALE 

ELM 

ELMGEN 

ETAGEN 

5JUOTNT 

REFGEN 

R0T3 

ROTGEN 

SCMOVE 

V2N0RM 

ZEROVC 

and  the  subroutine  ItCHPAR,  which  defines  machine-dependent  constants  (see  Section  11). 

The  subroutines  in  the  NPSOL  package  use  the  following  labelled  COMMON  areas: 

SOLMCH  (15  REAL  variables;  see  Section  11) 

S0L1CM  (3  INTEGER  variables) 

S0L3CM  (4  INTEGER  variables) 

S0L4CM  (10  REAL  variables) 

S0L1LP  (15  INTEGER  variables) 

S0L1NP  (30  INTEGER  variables) 

S0L2NP  (2  INTEGER  variables). 
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9.  DESCRIPTION  OF  THE  PRINTED  OUTPUT 

The  following  is  a  description  of  the  terse  line  printed  at  each  major  iteration  if  the  last  two 
digits  of  MSGLVL  >  5.  The  printout  from  the  QP  subroutines  is  described  in  Gill  et  a/.  (1983a). 
All  quantities  arc  evaluated  at  the  end  of  the  iteration. 


OBJECTIVE 


NCOLZ 


NORM  GFREE 


NORM  QTG 


NORM  ZTG 


CORD  H 


COND  T 


is  the  major  iteration  count,  k. 

is  the  number  of  minor  iterations  needed  to  solve  the  QP  subproblem. 

is  the  step  a*  taken  along  the  computed  search  direction. 

is  the  total  number  of  evaluations  of  the  problem  functions. 

is  the  value  of  the  objective  function,  /‘'(z*). 

is  the  number  of  bounds  in  the  predicted  active  set. 

is  the  number  of  linear  constraints  in  the  predicted  active  set. 

is  the  number  of  nonlinear  constraints  in  the  predicted  active  set. 

is  N  minus  the  number  of  constraints  in  the  predicted  active  set. 

is  the  norm  of  the  gradient  of  the  objective  function  with  respect  to  the 
free  variables  (not  printed  if  ORTHOG  is  .FALSE.). 

is  a  weighted  norm  of  the  gradient  of  the  objective  function  with  respect 
to  the  free  variables  (not  printed  ir  ORTHOG  is  .TRUE.). 

is  the  Euclidean  norm  of  the  projected  gradient. 

is  a  lower  bound  on  the  condition  number  of  the  Hessian  approximation, 
i.c.  a  bound  on  cond(//)  =  cond(/?r/£). 

is  a  lower  bound  on  the  condition  number  of  the  matrix  of  predicted 
active  constraints. 


NORM  C 


is  the  norm  of  the  vector  of  constraint  violations  and  residuals  of  the 
constraints  in  the  predicted  active  set. 


is  the  penalty  parameter  used  in  the  augmented  Lagrangian  merit  func- 
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CONY  is  a  four-letter  indication  of  the  status  of  the  four  convergence  tests; 

each  letter  is  “T”  if  the  test  is  satisfied,  and  “F”  otherwise.  The  four 
tests  indicate  whether:  (a)  the  projected  gradient  is  small;  (b)  the  active 
constraint  residuals  are  small;  (c)  the  multipliers  indicate  optimality;  (d) 
the  last  change  in  X  was  small. 

U  refers  to  the  quasi-Newton  update  of  R  to  obtain  a  new  estimate  of  the 

Hessian.  Uis  1  if  the  update  was  performed,  and  0  if  no  update  occurred. 

The  following  is  a  description  of  the  solution  output  of  NPSOL.  Note  that  names  are  automati¬ 
cally  assigned  to  each  variable  and  constraint. 

The  following  printout  is  given  for  each  variable  zy. 

VARIABLE  is  the  name  (VARBL)  and  index  j,  j  —  1  to  N,  of  the  variable. 

STATE  gives  the  state  of  the  variable  (FR  if  not  in  the  working  set,  EQ  if  in 

the  working  set  as  a  fixed  variable,  LL  if  in  the  working  set  at  its  lower 
bound,  and  UL  if  in  the  working  set  at  its  upper  bound).  If  VALUE  lies 
outside  the  upper  or  lower  bounds  by  more  than  FEATOL(j),  STATE  will 
be  or  respectively. 


VALUE 

LOVER  BOUND 

UPPER  BOUND 


is  the  value  of  the  variable  zy  at  the  final  iteration, 
is  the  lower  bound  BL (j)  specified  for  the  variable, 
is  the  upper  bound  B U(y)  specified  for  the  variable. 


LAGR  MULTIPLIER  is  the  value  of  the  Lagrange  multiplier  for  the  corresponding  bound 

constraint.  This  will  be  zero  if  STATE  is  FR.  If  X  is  optimal  and  STATE  is 
LL,  the  multiplier  should  be  non-negative;  if  STATE  is  UL,  the  multiplier 
should  be  non-positive. 

RESIDUAL  is  the  difference  between  the  variable  and  its  nearer  bound. 

The  following  printout  is  given  for  each  constraint. 

LINEAR  CONSTR  is  the  name  (LNCON)  and  index  t,  *  =  l  to  NCLIN,  of  a  linear  constraint. 


NONLMR  CONSTR 


is  the  name  (NLCON)  and  index  t,  i  =  1  to  NCNLN,  of  a  nonlinear  con- 
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STATE 


VALUE 


LOVER  BOUND 


UPPER  BOUND 


LAGR  MULTIPLIER 


RESIDUAL 


is  the  state  of  the  constraint  (FR  for  a  constraint  not  in  the  working  set, 
EQ  for  an  equality  in  the  working  set,  LL  for  an  inequality  constraint 
in  the  working  set  at  its  lower  bound,  UL  for  an  inequality  constraint 
in  the  working  set  at  its  upper  bound).  STATE  will  be  “++”  or  “ — ” 
respectively  if  VALUE  lies  outside  the  upper  or  lower  bounds  by  more 
than  its  feasibility  tolerance. 

is  the  value  of  the  constraint  at  the  final  point. 

is  the  specified  lower  bound  for  the  constraint. 

is  the  specified  upper  bound  for  the  constraint. 

is  the  value  of  the  Lagrange  multiplier.  This  will  be  zero  if  STATE  is  FR. 
If  X  is  optimal  and  STATE  is  LL,  the  multiplier  should  be  non-negative; 
if  STATE  is  UL,  the  multiplier  should  be  non-positive. 

is  the  residual  of  the  constraint  with  respect  to  its  nearer  bound,  i.e., 
the  difference  between  VALUE  and  the  nearer  of  the  two  bounds. 


10.  ERROR  RECOVERY 
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The  input  data  for  NPSOL  should  always  be  checked  (even  if  NPSOL  terminates  with  the 
value  INFORM  =  01).  Two  common  sources  of  error  are  uninitialised  variables  and  incorrect 
gradients,  which  may  cause  underflow  or  overflow  on  some  machines.  The  user  should  check  that 
all  components  of  A,  BL,  BU,  FEATOL  and  X  are  defined  on  entry  to  NPSOL,  and  that  OBJFUN  and 
CONFUN  compute  all  relevant  components  of  OBJGRD,  C  and  CJAC. 

The  present  version  of  NPSOL  contains  no  procedure  for  checking  the  computed  gradients. 
Incorrect  gradients  may  lead  to  termination  with  INFORM  =  2,  3  or  5. 

Other  error  conditions  may  arise  as  follows. 


Termination 


Recommended  Action 


Underflow 


Overflow 


INFORM  =  1 


If  the  machine  parameter  indicating  an  underflow  check  (WMACH(9))  is 
zero,  floating-point  underflow  may  occur  occasionally,  but  can  usually  be 
ignored.  To  avoid  underflow,  set  VMACH(9)  to  a  positive  value;  however, 
this  will  lead  to  a  noticeable  loss  of  efficiency.  If  underflow  continues  to 
occur  for  no  apparent  reason,  contact  the  authors  at  Stanford  University. 

If  the  printed  output  before  the  overflow  error  contains  a  warning  about 
serious  ill-conditioning  in  the  working  set  when  adding  the  j-th  con¬ 
straint,  it  may  be  possible  to  avoid  the  difficulty  by  increasing  the 
magnitude  of  FEATOL(j),  and  rerunning  the  program.  If  the  message 
recurs  even  after  this  change,  the  offending  linearly  dependent  constraint 
must  be  removed  from  the  problem.  If  overflow  occurs  in  one  of  the 
user-supplied  routines  (e.g.,  if  the  nonlinear  functions  involve  exponen¬ 
tials  or  singularities),  it  may  help  to  specify  tighter  bounds  for  some  of 
the  variables  (i.e.,  reduce  the  gap  between  appropriate  lj  and  uy).  If 
overflow  continues  to  occur  for  no  apparent  reason,  contact  the  authors 
at  Stanford  University. 

A  feasible  point  could  not  be  found  for  the  bounds  and  linear  constraints. 
This  exit  occurs  if  there  is  a  failure  in  the  LP  phase  of  any  QP  subproblem 
(see  Gill  et  af.,  1983a).  The  most  likely  reason  for  this  condition  is  that 
the  linear  constraints  and  bounds  arc  incompatible  or  inconsistent;  if 
so,  NPSOL  will  terminate  during  the  first  major  iteration.  In  order  for 
a  feasible  point  to  exist,  the  constraints  must  be  re-formulated,  or  the 
corresponding  components  of  FEATOL  must  be  re-defined,  as  discussed  in 
Gill  et  al.  (1983a).  Another  possibility  is  that  dependencies  among  the 
constraints  and  bounds  have  led  to  cycling  in  the  LP  phase;  this  will 
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always  be  the  case  if  NPSOL  terminates  with  INFORM  =  1  after  the  first 
major  iteration. 

A  sufficient  decrease  in  the  merit  function  could  not  be  attained  during 
the  final  line  search.  This  sometimes  occurs  because  an  overly  stringent 
accuracy  has  been  requested,  i.e.,  FTOL  is  too  small;  in  this  case  the 
final  solution  may  be  acceptable  despite  the  non-zero  value  of  INFORM 
(see  Gill,  Murray  and  Wright,  1981,  for  a  discussion  of  the  attainable 
accuracy).  If  the  projected  gradient  at  the  final  point  is  not  small,  the 
computed  gradients  may  be  incorrect.  Another  possibility  is  that  the 
search  direction  has  become  inaccurate  because  of  ill-conditioning  in  the 
Hessian  approximation  or  the  matrix  of  constraints  in  the  working  set; 
either  form  of  ill-conditioning  also  tends  to  be  reflected  in  large  values  of 
ITQP  (the  number  of  iterations  required  to  solve  each  QP  subproblein). 
If  the  condition  estimate  of  the  Hessian  (COND  H)  is  extremely  large,  it 
may  be  worthwhile  to  try  a  warm  start  at  the  final  point  with  COLD  set 
to  .FALSE.,  ISTATE  unaltered,  and  R  set  to  the  identity  matrix.  If  the 
matrix  of  constraints  in  the  working  set  is  ill-conditioned  (i.e.,  COND  T 
is  extremely  large),  it  may  be  helpful  to  run  NPSOL  with  relaxed  values 
of  the  components  of  FEATOL  corresponding  to  nearly  dependent  con¬ 
straints.  (Constraint  dependencies  are  often  indicated  by  wide  variations 
in  size  in  the  diagonal  elements  of  the  matrix  T,  whose  diagonals  will  be 
printed  if  the  last  two  digits  of  MSGLVL  >  30.) 

If  the  algorithm  appears  to  be  making  progress,  the  value  of  ITMAX  may 
be  too  small.  If  so,  increase  ITMAX  and  rerun  NPSOL  (possibly  using  the 
warm  start  facility).  If  the  algorithm  seems  to  be  “bogged  down”,  the 
user  should  check  for  incorrect  gradients  or  ill-conditioning  as  described 
above  under  INFORM  =  2.  Note  that  ill-conditioning  in  the  working  set  is 
sometimes  resolved  automatically  by  the  algorithm,  in  which  case  per¬ 
forming  additional  iterations  may  be  helpful.  However,  ill-conditioning 
in  the  Hessian  approximation  tends  to  persist  once  it  has  begun,  so  that 
allowing  additional  iterations  without  altering  R  is  usually  inadvisable.  If 
the  constraint  violations  have  not  been  significantly  reduced,  the  prob¬ 
lem  may  have  no  feasible  point. 

A  guaranteed  procedure  for  resolving  extremely  small  Lagrange  multi¬ 
pliers  has  not  been  included  in  NPSOL,  since  it  would  be  inherently  com¬ 
binatorial  (sec  Gill,  Murray  and  Wright,  1981,  for  further  discussion). 
In  some  cases,  the  difficulty  may  be  avoided  by  removing  certain  active 
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INFORM  =  5 


constraints  with  very  small  multipliers  from  the  problem,  and  rerunning 
NPSOL. 

With  exact  arithmetic,  the  search  direction  should  always  be  a  descent 
direction  for  the  merit  function.  If  this  value  of  INFORM  occurs,  the  com¬ 
puted  gradients  may  be  incorrect,  or  ill-conditioning  may  have  destroyed 
thf»  accuracy  of  the  search  direction.  The  user  should  check  for  these 
conditions  as  described  above  under  INFORM  =  2. 
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This  program  has  been  written  in  ANSI  (1966)  Fortran  and  tested  on  an  IBM  3081  computer 
using  the  WATFIV  Compiler,  Version  1,  Level  6.  All  subroutines  in  NPSOL  are  PFORT- 
compatible  (Ryder,  1974),  except  for  some  A2  Hollerith  specifications. 

At  the  beginning  of  NPSOL,  the  subprogram  MCHPAR  is  called  to  assign  various  machine- 
dependent  parameters.  These  parameters  are  stored  in  the  array  WMACH(15)  in  the  labelled  COMMON 
block  SOLMCH. 

The  specification  of  MCHPAR  is 

SUBROUTINE  MCHPAR 

REAL  VMACH 

COMMON  /SOLMCH/  VMACH  (15) 


The  first  eleven  components  of  the  REAL  array  WMACH  must  be  set  in  MCHPAR.  The  components 
of  VMACH  are  defined  as  follows. 


Definition 


VMACH(l) 

WMACH(2) 

VMACH(3) 

VMACH(4) 

VMACH(5) 

VMACH(6) 

VMACH(7) 

VMACH(8) 

VMACH(9) 


VMACH(IO) 

VMACH(ll) 


is  NBASE,  the  base  of  floating-point  arithmetic. 

is  NDIGIT,  the  number  of  NBASE  digits  of  precision. 

is  EPSMGH,  the  floating-point  precision. 

is  RTEPS,  the  square  root  of  EPSMCH. 

is  FLMIN,  the  smallest  positive  floating-point  number. 

is  RTMIN,  the  square  root  of  FLMIN. 

is  FLMAX,  the  largest  positive  floating-point  number. 

is  RTMAX,  the  square  root  of  FLMAX. 

is  UNDFLV,  which  specifies  whether  or  not  NPSOL  should  check  for 
underflow  in  certain  computations.  If  UNDFLV  =  0,  no  underflow 
checking  will  be  performed.  If  UNDFLV  is  set  to  a  positive  number, 
NPSOL  will  check  Tor  underflow  and  will  replace  too-small  quantities 
by  iero.  Note  that  NPSOL  will  run  faster  if  no  underflow  checking 
takes  place. 

is  NIN,  the  file  number  for  the  input  stream, 
is  NOUT,  the  file  number  for  the  output  stream. 
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The  following  version  of  MCHPAR  (which  is  provided  by  the  Systems  Optimization  Laboratory) 
contains  the  parameters  associated  with  double  precision  on  a  machine  in  the  IBM  370  scries. 
The  user  must  substitute  a  version  of MCHPAR  that  is  appropriate  1‘or  the  machine  to  be  used. 


SUBROUTINE  MCHPAR 
C 

DOUBLE  PRECISION  NMACH 
COMMON  /SOLHCH/  NMACH ( 15) 

C 

C  MCHPAR  MUST  DEFINE  THE  RELEVANT  MACHINE  PARAMETERS  AS  FOLLOWS. 

C  WMACH!  I  )  =  NBASE  =  BASE  OF  FLOATING-POINT  ARITHMETIC. 

C  NMACH ( 2 )  s  NDI6IT  =  NO.  OF  BASE  NMACH ( 1)  DI6ITS  OF  PRECISION. 

C  NMACH! 3)  =  EPSMCH  a  FLOATING-POINT  PRECISION. 

C  WMACH ( 4 1  a  RTEPS  a  SORT!  EPSMCH ) . 

C  NMACH! 5)  =  FLMIN  a  SMALLEST  POSITIVE  FLOATING-POINT  NUMBER. 

C  NMACH! 6)  =  RTMIN  a  SORT! FLMIN). 

C  NMACH! 7)  *  FLMAX  a  LARGEST  POSITIVE  FLOATING-POINT  NUMBER. 

C  NMACH! 8)  =  RTMAX  a  SORT! FLMAX). 

C  NMACH! 9)  a  UNDFLN  a  o.O  IF  UNDERFLOW  IS  NOT  FATAL.  *VE  OTHERWISE. 

C  NMACH! 10)  a  NIN  a  STANDARD  FILE  NUMBER  OF  THE  INPUT  STREAM. 

C  NMACH! 11)  a  NOUT  a  STANDARD  FILE  NUMBER  OF  THE  OUTPUT  STREAM. 

C 

INTEGER  NBASE.  NDIGIT .  NIN,  NOUT 

DOUBLE  PRECISION  DSQRT 
C 

NBASE 
NDIGIT 
WMACH! 1 1 
NMACH! 2) 

NMACH! 3) 

WMACH! 4) 

WMACH! 5) 

NMACH! 6) 

WMACH! 7) 

WMACH! 8) 

NMACH! 9) 

NIN 
NOUT 

WMACH! 10) 

NMACH! 1 1  ) 

C 

C -  IN  HATFIV,  ALLOW  UP  TO  100  UNDERFLOWS. 

C -  CALL  TRAPS  !  0.0,100  ) 

RETURN 
C 

C  END  OF  MCHPAR 
END 


=  16 
a  14 
a  NBASE 
a  NDIGIT 

a  WMACH! 1 )**( 1  -  NDIGIT) 
=  DSQRT! NMACH! 3 ) ) 
a  WMACH! 1 )««!  -62 ) 
a  DSQRT! NMACH! 5)) 
a  NMACH! 1 )*»61 
a  DSQRT! NMACH! 7 ) ) 

a  0.00*0 

a  5 
=  6 
a  NIN 
a  NOUT 


NPSOL/34 


11.  IMPLEMENTATION  INFORMATION 


The  values  of  NBASE,  NDIGIT,  EPSMCH,  FLMIN  and  FLMAX  for  several  machines  are  given  in  the 
following  table,  for  both  single  and  double  precision;  RTEPS,  RTMIN  and  RTMAX  may  be  computed 
using  Fortran  statements.  The  values  NIN  and  NOUT  depend  on  the  machine  installation. 

For  each  precision,  we  give  two  values  for  EPSMCH,  FLMIN  and  FLMAX.  The  first  value  is  a 
Fortran  decimal  approximation  of  the  exact  quantity;  use  of  this  value  in  MCHPAR  should  cause 
no  difficulty  except  in  extreme  circumstances.  The  second  value  is  the  exact  mathematical 
representation. 

Table  of  machine- dependent  parameters 


Variable 

IBM  360/370 
Single 

CDC  6000/7000 
Single 

DEC  10/20 
Single 

Univac  1100 

Single 

DEC  VAX 

Single 

| 

16 

2 

2 

2 

2 

NDIGIT 

6 

48 

27 

27 

24 

EPSMCH 

9.54E-7 

7.11E-16 

7.46E-9 

1.50E-8 

1.20E-7 

16-5 

2“« 

2-27 

2-26 

2-23 

FLMIN 

1.0E-78 

1.0E-293 

1.0E-38 

1.0E-38 

1.0E-38 

l6-«5 

2-975 

2-129 

2-m 

2-128 

Variable  IBM  360/370  CDC  6000/7000 


NDIGIT 


EPSMCH 


Double 


16 


14 


2.22D-16 

1613 


1.0D-78 

l8-«5 


Double 


2 


96 


2.53D-29 

2-»5 


1.0D-293 

2-975 


DEC  10/20 
Double 


2 


62 


2.17D-19 

2-62 


1.0D-38 

2-!29 


Univac  1100 
Double 


2 


61 


DEC  VAX 
Double 


2 


56 


2.78D-17 

2-55 


1.0D-38 

2-128 
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12.  EXAMPLE  PROGRAM  AND  OUTPUT 

This  section  contains  a  listing  and  the  computed  results  from  a  sample  main  program  that 
calls  NPSOL  to  solve  one  version  of  the  so-called  “hexagon”  problem  (a  different  formulation  is 
given  as  Problem  108  in  Hock  and  Schittkowski,  1081).  The  problem  is  to  determine  the  hexagon 
of  maximum  area  such  that  no  two  of  its  vertices  are  more  than  one  unit  apart  (the  solution  is 
not  a  regular  hexagon). 

All  constraint  types  are  included  (bounds,  linear,  nonlinear),  and  the  Hessian  of  the  Lagrangian 
function  is  not  positive  definite  at  the  solution.  The  problem  has  nine  variables,  non-infinite 
bounds  on  six  of  the  variables,  four  general  linear  constraints,  and  fifteen  nonlinear  constraints. 
The  objective  function  is 

F(x)  —  —  ZgZg  +  Xi*7  —  Xj*7  —  ZgZg  +  Z4Zg  +  ZgZg. 

The  bounds  on  the  variables  are 

xi  >0,  z6  >  0,  z«  >0,  Z7  >  0,  Zg  <  0,  and  z«  <  0. 

Thus, 

la  =  (  0,  -00,  ~oo,  —00,  0,  0,  0,  -00,  —  oo)r 
UB  =  (  +oo,  +00,  +00,  +00,  +00,  +00,  +00,  0,  0)r. 

The  general  linear  constraints  are 

zg  —  *1  ^  0,  Z3  —  zg  >  0,  Zg  —  Z4  >  0,  and  x4  —  x8  >  0. 


Hence, 


The  fifteen  nonlinear  constraint  functions  are 

ci(*)”*i+*s,  cg(z)-(z2-zi)a+(z7-zg)3,  c3(z)-(z3-z,)a+z2, 

c4(z)-(zi -Z4)a+(zg-Zg)2,  cg(z)  “  (zj  -  Zg)2  +  (zg  -  Xg)a,  Cg(z)-z|+z?, 

C7(*)-(*3  -*a)2  +*7,  Cg(z)-(z4  -Zg)a  +(z8  -Z7)3,  Cg(z)-(zg  -z5)a  +(z7  -  Zg)3, 

cio(*)“*3.  Cu(z)-(z4 -z3)a +zl,  Cig(z)-(z5-z3)a+z|, 

Ci3(*)-*a+*l,  cH(z)-(z4 -zg)2+(zg-zg)3,  C18(z)-zj+z5. 
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(For  most  applications,  it  would  be  preferable  to  replace  the  tenth  nonlinear  constraint  (z®  <  1) 
by  the  bounds  —  1  <  *3  <  1. 

The  nonlinear  constraints  are  all  of  the  form 

Ci{x)  <  1,  i  =  l,...,  15; 

hence,  all  components  of  lN  are  —00,  and  all  components  of  uN  are  1. 

The  starting  point  zo  is 

z0  =  ( .33,  .87,  1.1,  .67,  .33,  .33,  .67,  -.33,  -.67  )r, 

and  F(xo)  =  —1.4333  (to  five  figures).  The  optimal  solution  (to  five  figures)  is 

x  =  ( .060947,  .59765,  1.0,  .59765,  .060947,  .34377,  .5,  -.5,  -.34377  )T, 

and  F(x)  =  —1.34996.  (The  optimal  objective  function  is  unique,  but  is  achieved  for  other  values 
of  z.)  Six  nonlinear  constraints  are  active  at  x.  The  sample  solution  output  is  given  later  in  this 
section,  following  the  sample  main  program  and  problem  definition. 
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1 

2 

11 

15 


16 

17 

18 

19 

20 


21 


22 

23 

26 

25 


C 

C  EXAMPLE  moon  AM  FOR  SUBROUTINE  MP90L 
C  DOUBLE  PRECISION  VERSION  1.1.  APRIL  1083. 

C  THE  VALUES  OF  THE  PARAMETERS  EPSAF,  FTOL,  AND  FEATOL  ARE 
C  APPROPRIATE  FOR  A  MACHINE  WITH  A  PRECISION  OF  15  DECIMAL  DIBITS. 

INTEGER 
INTEGER 
INTEGER 
INTEGER 
INTEGER 

DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
OOUBLE  PRECISION 
LOGICAL 
EXTERNAL 
DATA 
* 

C 

C  SET  THE  DECLARED  ARRAY  DIMENSIONS. 

C  NROMA  >  THE  DECLARED  RON  DIMENSION  OP  A. 

C  NRONJ  =  THE  DECLARED  ROM  DIMENSION  OP  CJAC. 

C  NROMt  «  THE  DECLARED  RON  DIMENSION  OP  R. 

C  LIMORK  a  THE  LENGTH  OF  THE  INTEGER  MORK  ARRAY. 

C  LNORK  «  THE  LENGTH  OP  THE  OOUBLE  PRECISION  MORK  ARRAY. 

C 

NROMA  a  5 
MtONJ  a  28 
NROHR  a  10 
LIMORK  a  SO 
LNORK  a  1008 
C 

C  SEY  THE  APPROXIMATE  MACHINE  PRECISION. 

C 

EPSMCH  a  1.00-15 
C 

C  SET  THE  PROBLEM  DIMENSIONS. 

C  N  a  THE  NUMBER  OP  VARIABLES. 

C  NCLIN  a  THE  NUMBER  OP  GENERAL  LINEAR  CONSTRAINTS  (HAY  BE  01. 

C  NCNLN  a  THE  NUMBER  OP  NONLINEAR  CONSTRAINTS  (HAY  BE  0). 

C  NCTOTL  a  THE  TOTAL  NUMBER  OP  VARIABLES  AM)  CONSTRAINTS. 

C  (THE  ARRAYS  I3TATE,  BL,  BU.  CLAMBOA  MUST  BE  AT  LEAST 

C  THIS  LONG.) 

C 

N  »  8 
NCLIN  a  6 

uruiu  a  IS 

NCTOTL  «  N  ♦  NCLIN  ♦  NCNLN 
C 

C  ASSIGN  THE  DATA  ARRAYS. 

C  BOUNDS  .68.  BIGBND  MILL  BE  TREAYED  AS  PUIS  INFINITY. 

C  BOUNDS  .LI.  -  DIMM  MILL  RE  TREATED  AS  HMJB  INFINITY, 
r  NOUr  a  THE  UNIT  NUMBER  FOR  PRINTING. 


It  INFORM t  ITERt  ITMAX,  J.  LIMORK.  LNORK 
HSGLVLt  N»  NCLIN.  NCNLN,  NCTOTL 
NOUT.  MtOMA.  MtOMJ.  MtOMR.  NSTATE 
ISTATE( 28)  ' 

I WORK ( 50 ) 

BIGBND,  EPSAF,  EPSMCH,  RTEPS,  ETA,  FTOL,  OBJF 
A(5>9),  BL(28),  BU(28),  FEATOL( 28 ) 

020),  CJACI20.9),  CLAMDA( 28) 

OBJGRD( 9) ;  R(10,9>,  X(9) 

MDRK(IOOO) 

DSRRT 
ZERO,  ONE 

COLD,  FEAUN,  0RTH06 
OBJ FUN,  CONFUN 
ZERO  ,  ONE 
/0.00*0,  1.00*0/ 
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A  »  THE  SENERAL  CONSTRAINT  MATRIX. 

BL  >  THE  LONER  BOUNDS  ON  X.  A*X  AM)  C(X). 

BU  =  THE  UPPER  BOUNDS  ON  X.  A*X  AM)  CtX). 

X  *  THE  INITIAL  ESTIMATE  OF  THE  SOLUTION. 

NOUT  =  A 
BI68HD  *  1.00*10 
DO  30  J  =  1,  NCTOTL 
BL(  J)  =  -BIG8ND 
BUI J)  *  BI66ND 
30  CONTINUE 
BLII)  x  ZERO 
BLI5)  x  ZERO 
BLIA)  x  ZERO 
BLI7)  x  ZERO 

SET  LONER  BOUNDS  OF  ZERO  FOR  THE  FOUR  LINEAR  CONSTRAINTS. 

BL(IO)  x  ZERO 
BL1 11)  x  ZERO 
BL1 12)  x  ZERO 
BL1 13)  x  ZERO 

BUIS)  x  ZERO 
BUI 9)  x  ZERO 

SET  UPPER  BOUNDS  OF  ONE  FOR  ALL  IS  NONLINEAR  CONSTRAINTS. 


42 

DO  40  J  x  |4,  28 

43 

BUI  J)  >  ONE 

44 

C 

40  CONTINUE 

45 

XII)  x  .330*0 

44 

XI 2)  x  .470*0 

47 

XI 3)  x  1.10*0 

4B 

XI4)  x  .470*0 

49 

XI5)  x  .330*0 

50 

XI4)  «  .330*0 

51 

XI 7)  x  .470*0 

52 

XIO)  *  -.330*0 

53 

c 

XI 9)  x  -.470*0 

54 

00  40  J  >  1,  N 

55 

00  SO  I  X  1 ,  | 

54 

A( I,J)  x  Z1 

57 

50  CONTINUE 

5B 

40  CONTINUE 

59 

A( 1 ,1 )  x  -ONE 

40 

All, 2)  x  one 

41 

A(2,2)  x  -ONE 

42 

AI2.3)  *  ONE 

43 

A(3,3)  >  ONE 

44 

A(3»4)  >  -ONE 

45 

AI4.4)  x  ONE 

44 

A(4,5)  «  -ONE 

PRINT  THE  DATA. 

MUTE  (NOUT,  fitOO) 
00  70  I  x  1 ,  NCLIN 
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69  WRITE  INOUT >  22001  X,  IAII,J),  J=«  »H) 

70  70  CONTINUE 

71  WXTE  INOUT,  2300)  (BLtJ).  J*1 .NCTOTL) 

72  NRXTE  INOUT,  2400)  (BUIJ),  J=1 .NCTOTL) 

73  NRXTE  INOUT,  2500)  I  XIJ),  J=1»N> 

C 

c 

C  ALUM  UP  TO  50  MAJOR  ITERATIONS  TO  FIND  A  SOLUTION. 

C 

74  ITMAX  -  50 

C 

C  ASK  FOR  BRIEF  OUTPUT  EACH  MAJOR  ITERATION,  AM)  A  FULL  PRINT-OUT  OF 
C  THE  FINAL  SOLUTION. 

C 

75  MS6LVL  =  10 

C 

C  SET  THE  ABSOLUTE  PRECISION  OF  THE  OBJECTIVE  AT  THE  STARTING  POINT. 

C 

76  NSTATE  =  1 

77  CALL  OBJFUNt  2.  N,  X,  OBJF,  OBJGRD,  NSTATE  ) 

78  EPSAF  =  EPSMCH  «  DABS!  OBJF  ) 

C 

C  USE  A  SLACK  LINESEARCH. 

C  SET  THE  REQUIRED  NUMBER  OF  CORRECT  FIGURES  IN  THE  OPTIMAL  OBJECTIVE. 
C  THE  VALUE  CHOSEN  HERE  IFTOL  s  10  EPSMCH)  ASKS  FOR  ALMOST  FULL 
C  PRECISION  IN  OBJF. 

C 

79  ETA  s  0.90*0 

80  FTOL  s  10.00*0  *  EPSMCH 

C 

C  AT  THE  SOLUTION,  ANT  CONSTRAINT  HAY  BE  VIOLATED  BY  AS  MUCH  AS 
C  THE  SQUARE  ROOT  OF  THE  MACHINE  PRECISION. 

C 

81  RTEPS  s  DSQRTI  EPSMCH  ) 

62  DO  80  J  =  1 ,  NCTOTL 

83  FEATOL1 J)  3  RTEPS 

84  60  CONTINUE 
C 

C  A  COU)  START  IS  NEEDED  FOR  THE  FIRST  CALL  TO  NPSOL. 

C  START  THE  NONLINEAR  ITERATIONS  AT  A  POINT  THAT  IS  FEASIBLE  MXTH 
C  RESPECT  TO  THE  LINEAR  CONSTRAINTS  AN)  BOUNDS. 

C  USE  AN  ORTHOGONAL  FACTORIZATION  OF  THE  MATRIX  OF  CONSTRAINTS 
C  IN  THE  H0RKIN6  SET. 

C 

85  COLD  s  .TRUE. 

86  FEALIN  =  .TRUE. 

67  ORTHQG  *  .TRUE. 

C 

C 

C  SOLVE  THE  PROBLEM. 

C 

CALL  NPSOL!  ITMAX,  HSGLVL,  N, 

»  NCLIN,  NCNLN,  NCTOTL,  NROMA,  NRONJ,  NROMR, 

»  BXGBND,  EPSAF,  ETA,  FTOL, 

«  A,  8L,  BU,  FEATOL, 

*  CONFUN,  OBJ FUN,  COLD,  FEALIN,  ORTHOS, 

*  INFORM,  ITER,  ISTATE, 

«  C,  CJAC,  CLAMDA,  OBJF,  OBJGRD,  R,  X, 

*  WORK,  LXMORK,  HORK,  LMORK  ) 

r 


88 


Nl‘SOL/30 


12.  KXAMPMi  PROGKAM  AND  OUTPUT 


C  TEST  POT  AN  ERROT  CON)ZTXON. 

C 

69  XP  (XNPOTH  .OT.  6)  M  TR  966 

C 

C 

C  THE  FOLLOWING  IS  POT  ILLUSTRATIVE  PURPOSES  ONLY. 

C  ME  DO  A  MARK  START  NXTM  THE  FINAL  WORKING  SET  AND  R  OF  THE  PREVIOUS 
C  RUN,  BUT  WITH  A  SLIGHTLY  PERTURBED  STARTING  POINT. 

C 

90  DO  tOO  J  >  It  N 

91  XU)  a  XCJ)  ♦  0.050+0 

92  100  CONTINUE 
C 

C  RESET  THE  ABSOLUTE  PRECISION  OF  THE  OBJECTIVE  FUNCTION. 

C 

93  EPSAF  s  EPSNCH  *  DABS(  OBJF  ) 

C 

94  COLD  3  .FALSE. 

95  HS6LVL  3  5 

96  MRITE  (NOUT,  2600) 

97  MRITE  (NOUT,  2500)  (XU),  Jsl  ,N) 

C 

96  CALL  NPSOLt  XTMAX,  NS6LVL,  N, 

«  NCLIN,  NCNLN,  NCTOTL.  MNMA,  NRONJ,  MNMR, 

»  BIG8ND,  EPSAF.  ETA,  FTOL, 

«  A,  BL,  Ml.  FEATOL, 

*  CONFUN,  OBJ  FUN,  COLD.  FEALXN,  0RTN06, 

»  INFORM,  ITER,  1ST ATE, 

•  C.  CJAC,  C LAND A,  OBJF,  OBJSRD,  R,  X, 

*  I WORK,  LIMORK,  MORK,  LMORK  ) 

C 

99  IF  (INFORM  .ST.  0)  GO  TO  900 

100  STOP 
C 

C  ERROT  EXIT. 

C 

101  900  MRITE  (NOUT,  3000)  INFORM 

102  STOP 
C 

103  2100  FORMAT)/  12H  RONS  OF  A.) 

104  2200  FORMAT)/  (IX,  13,  4X,  9F6.2)) 

105  2300  FORMAT)/  14H  LONER  BOUNDS.  /  (IX,  1P7E10.2)) 

106  2400  FORMAT)/  14H  UPPER  BOUNDS.  /  (IX,  1P7E10.2)) 

107  2500  FORMAT(/  12H  INITIAL  X.  /  (IX,  7F10.2)) 

106  2600  FORMAT ( //46H  A  RUN  OF  THE  SAME  EXAMPLE  MXTH  A  HARM  START....) 

109  3000  FORMAT)/  32H  NPSOL  TERMINATED  M2TH  INFORM  3,  13) 

C 

C  END  OF  THE  EXAMPLE  PROGRAM  FOR  NPSOL. 

110  EM) 

111  SUBROUTINE  OBJFUNI  HOOE,  N,  X,  OBJF,  OBJGRO,  NSTATE  ) 

112  INTEGER  MODE,  N,  NSTATE 

113  OOUBLE  PRECISION  OBJF 

114  OOUBLE  PRECISION  X(N),  OBJGRD(N) 

C 

c  - — - 

C  OBJFUM  COMPUTES  THE  VALUE  AND  FIRST  DERIVATIVES  OF  THE  NONLINEAR 
C  OBJECTIVE  FUNCTION. 

C  - 

US  OBJF  «  -  X(t)«X(6)  ♦  X(1)«X(7)  -  X(3)«X(7)  -  X(5)*XI0) 

•  ♦  X(4MX(9)  *  X(3)«X(6> 
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116 

OBJOROin  * 

XI 7) 

117 

0BJGR0I2)  » 

“ 

X(6) 

IIS 

OBJGROt 3)  = 

- 

XI7) 

♦ 

XI8) 

119 

0BJ6RDI4)  = 

XI 9) 

120 

OBJGROIS)  * 

- 

XI 8) 

121 

0BJGR0I6)  = 

- 

XI 2) 

122 

0BJ6R0I7)  * 

- 

XI 3) 

♦ 

XII) 

123 

06JBRDI8)  a 

- 

XI5) 

♦ 

XI 3) 

124 

0BJ6R0I9)  « 

XI4) 

125 

RETURN 

c 

DO  OF  OBJFUN 

126 

END 

127 

SUBROUTINE  CONFUN! 

MODE.  NCNLN.  N, 

I2S 

INTEGER 

MODE,  NCNLN.  N. 

129 

DOUBLE  PRECISION 

XIN).  CINRQHJ). 

C 

c 

c 

c 

c 

c 

c 

c 


CONFUN  COMPUTES  THE  VALUES  AM)  FIRST  DERIVATIVES  OF  THE  NONUNEAR 
CONSTRAINTS. 

THE  ZERO  ELEMENTS  OF  JACOBIAN  MATRIX  ARE  SET  ONLY  ONCE.  THIS  OCCURS 
DURING  THE  FIRST  CALL  TO  CONFUN  INSTATE  =  t). 


130 

INTEGER 

I.  J 

131 

DOUBLE  PRECISION  ZERO.  TNO 

132 

OATA 

ZERO  .  TNO 

« 

SO. 00*0.  2.00*0/ 

133 

(p 

IF  INSTATE  . 

ME.  1 1  BO  TO  200 

134 

DO  120  J  3  1 

.  N 

135 

DO  110  I 

3  1,  NCNLN 

136 

CJACII.J)  3  zero 

137 

110  CONTINUE 

138 

120  CONTINUE 

139 

V 

200  Cl  1 ) 

X1 1  )»*2  *  XI6HW2 

140 

CJACI1.1) 

TN03XI1) 

141 

CJACtt.6) 

TN0«Xt6) 

142 

CI2) 

1X12)  -  XI1))M2  * 

143 

CJACI2.1 ) 

-  TN03IXI2)  -  XII)) 

144 

CJACt 2.2) 

TNWMXI2)  -  XII )) 

145 

CJACI2.6) 

-  TNO«IXI7)  -  XI6) ) 

146 

CJACI2.7) 

TN03IXI7)  -  XI 6 ) ) 

147 

w 

CI3) 

IXC3)  -  XII) HHI*  ♦ 

148 

CJACt 3.1) 

-  TM03IXI3)  -  XIII) 

149 

CJACI3.3) 

TN0«IX(3)  -  XII )) 

150 

CJACt  3.6) 

TM03XI6) 

151 

w 

CIA) 

1X11)  -  XI4))««2  * 

152 

CJACI4.1 ) 

THO«IXI1 )  -  XI 41) 

153 

CJACI4.4) 

-  TNO*IX1 1 )  -  XI4)) 

154 

CJACI4.6) 

TN03IXI6I  •  XI8)) 

155 

CJACt 4.81 

-  TN03IXI6)  -  XIOI) 

1X17)  -  XI6>)»*2 


X(6)«*2 


1X16)  -  X(8>)**2 


r 
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154 

CIS) 

8 

157 

CJACIS.1) 

B 

158 

CJACI5.5) 

X 

159 

CJACI 5,6) 

8 

160 

C 

CJACI 5, 9) 

B 

161 

CIS) 

B 

162 

CJACI  6, 2) 

B 

163 

C 

CJACI 6, 7) 

3 

164 

CI7) 

B 

165 

CJACI 7,2 ) 

8 

166 

CJACI 7,3) 

8 

167 

c 

CJACI 7,7) 

8 

168 

CIS) 

8 

169 

CJACI8.2) 

8 

178 

CJACI 8.4) 

3 

171 

CJACI 8,7) 

8 

172 

c 

CJACI 8,8) 

8 

173 

C19) 

8 

174 

CJACI 9.2) 

8 

175 

CJACI9.5) 

8 

176 

CJACI 9,7) 

a 

177 

c 

CJACI 9,9) 

8 

178 

CI10) 

X 

179 

c 

CJACI 10,3) 

8 

180 

CI11) 

X 

181 

CJAC1 11 ,3) 

8 

182 

CJACI 11. 4) 

8 

183 

c 

CJACI11.8) 

8 

184 

CMC) 

S 

185 

CJACI 12,3) 

S 

186 

CJACI 12.5) 

S 

187 

c 

CJACI 12,9) 

X 

188 

Cl  13) 

S 

189 

CJACI 13.4) 

B 

190 

c 

CJACI 13,8) 

8 

191 

Cl  14) 

S 

192 

CJACI 14,4) 

8 

193 

CJACI 14,5) 

B 

194 

CJACI 14,8) 

B 

195 

c 

CJACI 14,9) 

S 

196 

CI15) 

S 

197 

CJACI 15,5) 

S 

198 

CJACI 15,9) 

8 

1X11)  -  XIS) )*«2  ♦  (X(6)  - 

1X0*1X11)  -  XI5)) 

1XD*IXI  1 )  -  XISII 
1X0*1X16)  -  X19) ) 

1X0*1X16)  -  X( 9) J 

XI2>«*2  ♦  XI7)**2 
1XD*XI2> 

1X0*XI7) 

(X( 3)  -  XI2)>**2  ♦  XI7)**2 

ixo*ixt3>  -  xitn 
1X0*1X13)  -  Xltll 
TX0*XI7) 

ixi4)  -  xmw  ♦  txts)  - 

1X0*1X14)  -  XI 2 ) ) 

1X0*1X14)  -  XI2)) 

1X0*1X18)  -  XI 7)) 

1X0*1X18)  -  XI 7) ) 

1X12)  -  X(5)MH*2  *  1X17)  - 

1X0*1X12)  -  XIS)) 

1X0*1X12)  -  XfS)) 

1X0*1X17)  -  XI 9) ) 

1X0*1 XI7)  -  XI9)) 

X(3)M2 

1X0*XI3) 

1X14)  -  XI3)>**2  ♦  XI8)**2 

1X0*1X14)  -  XI 3)) 

1X0*1X14)  -  XI3)) 

TKD*XI8) 

1X15)  -  XI 3) )**2  ♦  XI9HW2 

1X0*1X15)  -  XI 3) ) 

1X0*1X15)  -  XIS) ) 

TM0*X19) 

XI 4  )**2  ♦  XI8)**2 
1M0*XI4) 

1M0*XI8) 

IXI4)  -  XI5))**2  ♦  1X19)  - 
7X0*1X14)  -  XI 5)) 

1X0*1X14)  -  XIS)) 

1X0*1X19)  -  XI8)) 

1X0*1X19)  -  XIS) ) 

XI5)**2  *  XI9)**2 
1X0«XI5) 

1X0*X(9) 


XI 9) )**2 


X17))**2 


X19))**2 


XIS)  )**2 
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ROMS  OF  A. 


1  -1 .00 

1. 

.00  0.00 

0.00 

2 

0.00 

-1. 

.00  1.00 

0.00 

3 

0.00 

0. 

.00  1.00 

-1.00 

4 

0.00 

0. 

.00  0.00 

1.00 

LONER  BOUNDS. 

0.000-01 

-1.000 

10 

-1.000  10  -1.000  10 

-1.000  10 

-1 .000 

10 

0.000-01 

0.000-01 

-1.000  10 

-1 .000 

10 

-1.000  10  • 

-1.000  10 

-1.000  10 

-1.000 

10 

-1.000  10  - 

-1.000  10 

UPPER  BOUNDS. 

1 .000  10 

1.000 

10 

1.000  10 

1.000  10 

0.000-01 

0.000- 

-01 

1.000  10 

1.000  to 

1.000  00 

1.000 

00 

1.000  00 

1.000  00 

1 .000  00 

1.000 

00 

1.000  00 

1.000  00 

INITIAL  X. 

0.33 

0.67 

1.10 

0.67 

-0.33  -0.67 


0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

-1.00 

0.00 

0.00 

0.00 

0.00 

0, 

.000- 

-01 

0.000- 

-01 

0. 

.000- 

-01 

0 

.000-01 

0.000- 

-01 

-1, 

.000 

10 

•1 

.000 

10 

-1.000 

10 

-1, 

.000 

to 

-1 

.000 

10 

-1.000 

10 

-1, 

.000 

10 

1 

.000 

10 

1.000 

10 

1, 

.000 

to 

1 

.000 

10 

1.000 

10 

1 

.000 

00 

1 

.000 

00 

1.000 

00 

1 

.000 

00 

1 

.000 

00 

1.000 

00 

1 

.000 

00 

0. 

.33 

0. 

.33 

0. 

.67 

WORKSPACE  PROVIDED  IS  INI  501.  H<  tOOO). 

TO  SOLVE  PROBLEM  HE  NEED  IH(  18).  HI  800). 


ITN  I TW 

STEP  HUHF  OBJECTIVE  BM) 

LC 

NC  NCOLZ  NORN  6FREE 

NORM  ZTB 

CCND  H 

COND  T 

NORN  C 

RHO  CONV  U 

0 

— 

0.00-01 

1  -1.44000  00 

— 

— 

-- 

— 

2.050 

00 

— 

-- 

— 

9.360-01 

— 

— 

- 

1 

a 

4.00-01 

3  -1 .47230  00 

0 

0 

5 

4 

2.080 

00 

9.770-02 

1.60 

00 

3.30  00 

6.980-01 

0.00-01 

FFFF 

2 

8 

1.00 

00 

4  -1.34230  00 

0 

0 

4 

5 

2.000 

00 

1 .670-01 

2.20 

00 

1.50 

00 

6.910-02 

0.00-01 

FFTF 

3 

1 

1.00 

00 

5  -1.33970  00 

0 

0 

4 

5 

2.050 

00 

1 .730-01 

6.90 

00 

1.40 

00 

4.570-02 

0.00-01 

FFTF 

4 

3 

1.40 

00 

7  -1.35470  00 

0 

0 

6 

3 

2.090 

00 

1.280-01 

6.90 

00 

3.00 

00 

1.110-01 

0.00-01 

FFTF 

0 

5 

1 

1.00 

00 

8  -1.35600  00 

0 

0 

6 

3 

2.060 

00 

3.210-02 

5.50 

00 

3.00 

00 

1 .690-02 

0.00-01 

FFTF 

6 

1 

1.00 

00 

9  -1.34980  00 

0 

0 

6 

3 

2.050 

00 

1 .080-02 

5.60 

00 

3.00 

00 

2.370-04 

0.00-01 

FFTF 

7 

1 

1.00 

00 

10  -1.34990  00 

0 

0 

6 

3 

2.050  00 

7.690-03 

7.50 

00 

3.00 

00 

1.290-04 

0.00-01 

FFTF 

8 

1 

1.00 

00 

11  -1.35000  00 

0 

0 

6 

3 

2.050 

00 

6.060-03 

1.50 

01 

3.00 

00 

2.430-04 

0.00-01 

FFTF 

9 

1 

1.00 

00 

12  -1.35000  00 

0 

0 

6 

3 

2.050 

00 

2.370-03 

2.80 

01 

3.00 

00 

1 .050-04 

0.00-01 

FFTF 

10 

1 

1.00 

00 

13  -1.35000  00 

0 

0 

6 

3 

2.050 

00 

4.150-04 

3.00 

01 

3.00 

00 

4.850-06 

0.00-01 

FFTF 

11 

1 

1.00 

00 

14  -1.35000  00 

0 

0 

6 

3 

2.050 

00 

4.570-05 

2.40 

01 

3.00 

00 

8.110-08 

0.00-01 

FFTF 

12 

1 

1.00 

00 

15  -1.35000  00 

0 

0 

6 

3 

2.050 

00 

7.400-06 

2.30 

01 

3.00 

00 

3.170-09 

0.00-01 

FTTF 

13 

1 

1.00 

00 

16  -1.35000  00 

0 

0 

6 

3 

2.050 

00 

7.430-07 

2.30 

01 

3.00 

00 

3.400-11 

0.00-01 

FTTF 

14 

1 

1.00 

00 

17  -1.35000  00 

0 

0 

6 

3 

2.050 

00 

2.890-08 

1.80 

01 

3.00 

00 

6.070-13 

0.00-01 

TTTF 

15 

1 

1.00 

00 

18  -1.35000  00 

0 

0 

6 

3 

2.050 

00 

1 .810-09 

1.40 

01 

3.00 

00 

9.440-15 

0.00-01 

TTTT 

EXIT  HP  PHASE.  INFORM  *  0  HAJXTS  ■  15  NFEVAL  «  18  NCEVAL  ■  18 


VARIABLE 

STATE 

VALUE 

LONER  BOUND 

VARBL 

1 

FR 

0.60946650-01 

0.0000000 

VARBL 

2 

FR 

0.5976493 

NONE 

VARBL 

3 

FR 

1.000000 

NONE 

VARBL 

4 

FR 

0.5976493 

NOME 

VARBL 

5 

FR 

0.60946650-01 

0.0800000 

UPPER  BOUM)  LA8R  MULTIPLIER  RESIDUAL 

NONE  0.0000000  0.60950-01 

NONE  8.0000000  0.10000  11 

0.0000000  0.1000D  II 

0.0000000  0.10000  11 

O.cmooooo  0.60950-01 


NONE 
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VARBL 

6 

FR 

0.3937715 

8.0000000 

NONE 

0.0000000 

0.3936 

VARBL 

7 

FR 

0.5000000 

0.0000000 

NONE 

0.0000000 

0.5000 

VARBL 

B 

FR 

-0.5000000 

NONE 

0.0000000 

0.0000000 

0.5000 

VARBL 

9 

FR 

-0.3937715 

NONE 

0.0000000 

0.3938 

LINEAR  CONSTR 

STATE 

VALUE 

LONER  BOUND 

UPPER  BOUM) 

LA6R  MULTIPLIER 

RESIDUAL 

INC  ON 

1 

FR 

0.5397027 

0.0000000 

NONE 

0.0000000 

0.5367 

INC  ON 

2 

FR 

0.9023507 

0.0000000 

NONE 

0.0000000 

0.9029 

LNCflN 

3 

FR 

0.9023507 

0.0000000 

NONE 

0.0000000 

0.9029 

LNCON 

9 

FR 

0.5367026 

0.0000000 

NONE 

0.0000000 

0.5367 

NONLIOI  CONSTR 

STATE 

VALUE 

LONER  BOUND 

UPPER  BOUND 

LAGR  MULTIPLIER 

RESIDUAL 

NLCON 

1 

FR 

0.1218933 

NONE 

1.000000 

0.0000000 

0.6781 

NLCON 

2 

FR 

0.3129571 

NONE 

1.000000 

0.0000000 

0.6875 

NLCON 

3 

UL 

1 .000000 

NONE 

1.000000 

-0.6318906D-01 

-0.92190-19 

NLCON 

9 

UL 

1.000000 

NONE 

1.000000 

-0.3202625 

-0.3997D-19 

NIC ON 

5 

FR 

0.9727152 

NONE 

1.000000 

0.0000000 

0.5273 

NLCON 

6 

FR 

0.6071697 

NONE 

1.000000 

0.0000000 

0.3926 

NLCON 

7 

FR 

0.9118861 

NONE 

1.000000 

0.5881 

NLCON 

a 

UL 

1.000000 

NONE 

1.000000 

-0.1992983 

-0.3997D-19 

NLCON 

9 

UL 

1.000000 

NONE 

1 .000000 

-0.3202625 

-0. 99910-19 

NLCON 

10 

UL 

1.000000 

NONE 

1.000000 

-0.3937715 

0.0000 

NLCON  11 

FR 

0.9118861 

Nora 

1.000000 

0.0000000 

0.5881 

NLCON 

12 

UL 

1.000000 
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This  report  forms  the  user's  guide  for  version  1.1  of  SOL/NPSOL,  a  set  of 
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NPSOL  uses  a  sequential  quadratic  programming  (SQP)  algorithm,  in  which  the 
search  direction;  is  the  solution  of  a  quadratic  programming  (QP)  sub¬ 
problem.  The  algorithm  treats  bounds,  linear  constraints  and  nonlinear 
constraints  separately.  The  Hessian  of  each  QP  subproblem  is  a  positive- 
definite  quasi-Newton  approximation  to  the  Hessian  of  an  augmented 
Lagrangian  function.  The  steplength  at  each  iteration  is  required  to 
produce  a  sufficient  decrease  in  an  augmented  Lagrangian  merit  function. 
Each  QP  subproblem  is  solved  using  a  quadratic  programming  package  with 
several  features  that  improve  the  efficiency  of  an  SQP  algorithm. 
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